function [k logk] = computeNormalisationConstant(S, tgt, tol)

if nargin < 3
    tol = 1e-6;
    if nargin < 2
        tgt = 1;
    end
end

% clamp distances to above 0
S = max(S, tol);

% initialise search range
min_logk = 0;
while computeZ(min_logk, S) <= tgt
    min_logk = min_logk-1;
end
max_logk = min_logk;
while computeZ(max_logk, S) >= tgt
    max_logk = max_logk+1;
end
min_logk_Z = computeZ(min_logk, S);
max_logk_Z = computeZ(max_logk, S);
assert(min_logk_Z >= tgt);
assert(max_logk_Z <= tgt);

% estimate log(k) by binary search
count = 0;
logk = 0.5*(min_logk+max_logk);
logk_Z = computeZ(logk, S);
while abs(logk_Z - tgt) > tgt*tol
    count = count + 1;
    if logk_Z < tgt
        max_logk = logk;
    else
        min_logk = logk;
    end
    logk = 0.5*(min_logk+max_logk);
    logk_Z = computeZ(logk, S);
end

k = exp(logk);



function R = computeRD(logk, SSD)

R = exp(-exp(logk)*SSD);



function Z = computeZ(logk, SSD)

R = computeRD(logk, SSD);

Z = sum(R(:));
